# Set parameters for abundance distrilim_unif <-100par1 <-c(0,log(10), log(3), 0) # the median is exp(par1)par2 <-c(lim_unif,log(1.5), log(1.5), log(10)) # CV is more or less sqrt(exp(sigma^2)-1)names(par1) <- abund_distrinames(par2) <- abund_distri
Below are the abundance distribution functions:
Code
x <-seq(0, lim_unif, by =0.01)glist <-list()for (t innames(par1)) {if (t =="1_unif") { y <-dunif(x, par1[t], par2[t]) tit <-paste0(t, " (min = ", par1[t], ", max = ", par2[t], ")") } else { tit <-paste0(t, " (meanlog = ", round(par1[t], 2), ", sdlog = ", round(par2[t], 2), ")") y <-dlnorm(x, meanlog = par1[t], sdlog = par2[t]) } df <-data.frame(x = x, y = y) g <-ggplot(df) +geom_line(aes(x = x, y = y)) +ggtitle(tit) +theme_linedraw() glist[[t]] <- g}wrap_plots(glist)
# Select the experiments for which delta should varyexp_vary_delta <-"ninter"
Code
# Create a dataframe of legends for the plotslegend_df <-data.frame(exp =unique(paramdf$exp))legend_df$legend <-c("Number of interactions","Species abundance distribution","Mean consumers niche breadth","Standard deviation of consumers niche breadth")
Simulations
perform_simu is FALSE so the simulations will be loaded from the files simres_delta0.rds, simres_delta0.2.rds, simres_delta1.rds.
Code
simres_list <-list()for (i in1:length(delta_list)) { delta <- delta_list[i]print(paste("Delta", delta, "----------------")) delta_name <- delta_names[i]if (perform_simu) { simres <-list("cor"=list("mean_cor"=data.frame(),"sd_cor"=data.frame(),"sd_area_cor"=data.frame()),"example"=list()) rtime <-system.time({if (delta == delta_def) {# Do all simulations jchoice <-1:nrow(paramdf) } else {# Do only the ones in delta_vary jchoice <-which(paramdf$exp %in% exp_vary_delta) }for(j in jchoice) {# Simulation parameters --- parms <- paramdf[j, ]print(paste0("Experiment ", parms$exp, " (", parms$id, ") ---"))for(repi in1:nrep) {print(paste("Rep", repi))# parms$rep <- repi # For debugging# Simulate data ---# Generate consumer and resource abundances # These abundances must be unsorted, because in the simulation# algorithm we sort the resource following their traits and if the abundances# are sorted as well it creates a correlationif (parms$abund_distri =="1_unif") { consumerab <-runif(parms$nconsumer, par1[parms$abund_distri], par2[parms$abund_distri]) resourceab <-runif(parms$nresource, par1[parms$abund_distri], par2[parms$abund_distri]) } else { consumerab <-rlnorm(parms$nconsumer, meanlog = par1[parms$abund_distri], sdlog = par2[parms$abund_distri]) resourceab <-rlnorm(parms$nresource, meanlog = par1[parms$abund_distri], sdlog = par2[parms$abund_distri]) }# Simulate data sim <-compas2d(nsp = parms$nconsumer, nsite = parms$nresource,le.grad = le.grad, col_prefix ="b",row_prefix ="p",rowvar_prefix ="tp",remove_zeroes =TRUE,col_abund = consumerab,row_abund = resourceab,ninter = parms$ninter,ratio.grad = ratio.grad,mean.tol = parms$mean_tol, sd.tol = parms$sd_tol, delta = delta,buffer = buffer)# Format simulated data --- comm <-data.frame(sim$comm) consumer_niche <- sim$colvar resource_traits <-data.frame(sim$rowvar)# Multivariate analyses ---# CA ca <-dudi.coa(comm, scannf =FALSE, nf =min(dim(comm)))# Reciprocal scaling rec <-reciprocal.coa(ca)# Measure correlations between mean positions --- res_niches <-get_niches(rec, consumer_niche = consumer_niche, resource_traits = resource_traits,comm = comm,rowname ="resources", colname ="consumers") res <-get_cor(res_niches) ressim <-add_params(res, parms)# Get eigenvalues neig <-10# We fix the number of eigevalues to 10, in case there are less eigenvalues there are NAs in the columns ca_eig <-as.data.frame(t(ca$eig[1:neig]))colnames(ca_eig) <-paste0("l", 1:neig) ca_eig[, names(parms)] <- parms# Chi-squared test sum_eig <-sum(ca$eig, na.rm =TRUE) n_comm <-sum(comm) testval <-chisq.test(comm, simulate.p.value =TRUE) pval <- testval$p.value# Add these to the results ca_eig$sum_eig <- sum_eig ca_eig$n_comm <- n_comm ca_eig$pval <- pval# Store results simres$cor <-combine_cor(simres$cor, ressim) simres$eig <-rbind(simres$eig, ca_eig) }# Keep one set of simulated data simres$example[[parms$id]] <- sim } })# Save resultsaveRDS(simres, file = simu_files[[delta_name]])# Store result in list simres_list[[delta_name]] <- simresmessage(paste0("Output of system.time for delta = ", delta, ":"))print(rtime)message("The simulation took ", round(rtime[3]/60, 2)," minutes to run.") } else { # Use stored objectfor(delta in delta_list) { simres_list[[delta_name]] <-readRDS(simu_files[[delta_name]]) } }}